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Abstract 

It is estimated that large quantities of lunar planetary pictures are 
to he obtained from forthcoming experiments. The processing of pictures through 
digital computer techniques offers one possibility for handling the large number 
of pictures which are to result from these experiments. This report is the final 
report on a theoretical and experimental study of techniques for processing 
pictures . 

An outline is presented of possible linear and non-linear theoretical 
work which can be made applicable to processing lunar and planetary pictures. 

The theoretical material of this report is limited to linear processing. The 
material deals with exploiting the geometry of pictures principally through the 
assumption of statistically stationary properties under translation and rotation. 
Pictures are defined as functions of one index where the index is a 2 -dimensional 
vector corresponding to the 2-dimensional coordinates of the picture. This nota- 
tion makes it possible to handle pictures in the conventional 1-dimensional nota- 
tion of classical processing theory while preserving the 2-dimensional properties 
of the picture. The report deals with representation of the pictures, linear 
processing, matrix products, correlation functions, stationary and symmetry pro- 
perties optimum linear processing, trivial processing, orthogonal preprocessing, 
minimum error, minimization with a constraint, preprocessing with prediction, 
quantization, and convergence of iterative computations. 

The experimental work of this report tests the theoretical results for 
complexity, usefulness, and correctness. An experimental capability was developed 
capable of producing pictures with up to 750 elements in 6 shades of brightness 
corresponding to the 6 faces of the cube onto which the shades of brightness were 
pasted. This capability was then used in the experimental representation of the 
Craters Eratosthenes and Archimedes in 15 x 19 fields of elements. Artificial 
pictures were also constructed having controlled correlation functions. A number 
of experiments were performed with these pictures to extract in an optimum manner 
another desired picture . A variety of FORTRAN programs are included for eventually 
performing the processing computations on a digital computer. 


1 
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Introduction 

The successful extraction of significant data from lunar and planetary 
television pictures requires the survey and comprehension of presently known tech- 
niques as well as the development of new techniques particularly directed at tele- 
vision processing. This report deals with theoretical and experimental work in 
the field of processing pictures. Much of the material of the report comes from 
3 previous reports'*"^. 

The theoretical work is directed at producing techniques for extracting 
data either before transmission or subsequent to transmission# Pre-processing of 
television pictures has the possibility of minimizing the transmission channel re- 
quirements or maximizing the use of the transmission channel when the channel is 
fixed and inflexible. Subsequent processing of television pictures after trans- 
mission permits an organization of the data into a form which is more meaningful 
as a final product and also permits hypothesis testing for the extraction of hy- 
pothetical data from the pictures. 

The quantity of data to be gathered from lunar picture experiments re- 
quires both automatic reduction and the automatic design of prototype reduction 
programs. Many . of the pictures will be used solely for mapping areas where only 
rough maps are presently available . Scientific efforts, however, will be directed 
at measuring the density distribution of crater radii and typing the debris around 
the craters as well as within the crater. Knowledge of the crater structure has 
the possibility of distinguishing between comet and asteroid impacts as well as 
the possibility of the recognition of volcanic craters. Similar crater studies 
will be performed on the scientific experiments directed at the Planet Mercury. 
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Geological studies will "be interested in the uniformity and roughness of the 
surface structure as well as the classification of the surface geology. 

Theoretical Studies 

Television pictures are normally thought of as a sequence of pictures 
capable of displaying moving objects. The use of television pictures in lunar 
and planetary experiments is much broader than this concept. Many of the pic- 
tures will be individually distinctive and will require appropriate processing. 
Some of the pictures will be related to others through observation of substan- 
tially the same experiment from different directions and different times . These 
sets of pictures require simultaneous processing which is considerably different 
than that required for the processing of a sequence of pictures of moving objects. 

Many optimum processing techniques are presently known for the proc- 
essing of telemetry. Very few of these methods have ever been used in the proc- 
essing of pictorial data. The following theoretical areas develop a number of 
techniques directly applicable to processing television pictures with particular 
emphasis on the use of. digital computation facilities for processing pictures 
both before and after transmission. The areas are listed in order of increasing 
difficulty. ■ 

Linear Special Filtering based on previously measured correlation 
functions can be used to preprocess data to reduce its quantity as well as to post 
process data to remove random irregularities such as noise and to enhance various 
effects for possible later detection. Most of the theoretical work presented in 
this report deals with this area. Some of the original work was done by Bell Tele 
phone Laboratories and is reported in references 7, 13, 14 and 24. Much of the 
theoretical work in 1 dimension as presented in references 11 and 1 6 is changed 
into 2 dimensions in this report. The concepts of rotation and translation have 
an expanded meaning in 2 dimensions. 
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Linear Learning Processes introduce the concept of measuring para- 
meters from sets of hypothetically characteristic data. The development of an 
intelligent learning process has the possibility of replacing the heuristic 
guessing which has normally accompanied picture processing in the past. The 
ideas of References 1, and k have the possibility of being applicable. Other 
approaches are contained in References 8, 15, and 25. 

Nonlinear Spat ial Fil tering and Learning is directed at processing the 
inherent nonlinear properties of a picture. The shades -of -gray scale is defined 
as the logarlthum of the brightness scale to correspond to the human preception 
of relative changes in the grey levels rhther than the absolute changes in the 
levels. The brightness scale also has an upper and lower bound. In order to 
handle these nonlinear effects and possibly others, it is necessary to seriously 
study techniques of nonlinear processing. Much of the basic theory is contained 
in Reference 27- Examples of the implementation of this theory will be found in 
References, 3, 6, 9, 17, and 20. 

Sequential Filtering of Pictures which are correlated with one another 
in sequence is an expansion of the amount of data used in the processing of a 
single point of a picture. The theory is applicable to scanning systems where 
it may be desirable to process the movement and change of objects in the pictures. 
The material of References 12 and 23 has possible application. 

Extraction of the presence or absence of data in a picture from only a 
hypothetical knowledge of what the data is to look like is an area particularly 
directed toward the discrete recognition or rejection of the presence of data in 
the picture under observation. The techniques of this area are also applicable 
to control of experiments and the acceptance of valid commands. The material of 
references 5, 8, and 18 can be applied to this subject. 
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Linear Analog Techniques provide an alternate means of processing pic- 
tures. Most theoretical work results in discrete processing on a digital com- 
puter. A number of multi -dimensional analog techniques also exist which have 
the possibility of producing simple linear equipment which may work at consider- 
ably higher speeds than can be obtained with digital equipment. The original 
work is contained in Reference 26. Further work is in References 2 , and 28. 

The following theoretical and experimental work develops and tests 
a number of linear techniques to handle the two-dimensional aspects of pictures 
so that the goemetry of the picture can be retained in two dimensions rather in 
one dimension where a considerable amount of theory is presently well known. 

Representation of a Picture 

The following theoretical development has 2 objectives: 1. to rep- 

resents picture and its processing in classical vector and matrix notation, 
and 2. to preserve the 2-dimensional properties of the picture. 

A sampled picture is represented by the row vector x having elements 
x^,^ corresponding to the samples of the picture where the indicies 1 = (l,l) 
and i = (a,b) are 2 -dimensional vectors. In classical notation, it is customary 
to write the elements of a row vector sequentially in one row. However, in the 
notation of this report a row vector is written in a geometrical form, each of 
the elements of the index vector indicating the geometrical position of the cor- 
responding row vector element . The elements of the representation of x then 
correspond geometrically to the picture that they represent. A possible rep- 


resentation of a row vector is then 

x * M = K] ’ 

\l 

X 12 

X 13 

x l4 


X 21 

X 22 

X 23 

x 24 


x 3l 

X 32 

X 33 

X 34 


Vfl 
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Linear Processing 

A second picture y can be obtained from the picture x by linearly 
weighting and adding the elements of x to form the elements of the picture y. 
Notationally, the new picture is computed by the matrix formula 


y = xh + k 


where 



h. . + 



The summation is over all the values of the index i. In reality, this is a 
double summation. It is represented here as a single sum because of its simpli 
city and its similarity to the classical matrix notation of the matrix product. 


Linear processing can be used to produce an isolated point in a pic- 
ture or to approximate some desired result based on the data contained in the 
processed picture. 


Matrix Transposition 


The row vector x is a degenerate form of a matrix h having elements 


h. . where the indices 1 and j are 2 -dimensional vectors. 
1 > 3 


The transpose h of 


the matrix h . has elements h 


3,1 


L h i,j] 


h l 


= [h. .1 


The transpose x^ of the row vector x = £x^ J is defined to be a column vector. 
Note that the indices have been transposed, not the elements of the vector rep- 
resenting the indices. The representation of the column vector x^ is the same 
as for the row vector x. 


X 11 

C\J 
i — i 
* 

X 13 

x ll 

X 21 

X 22 

X 23 

X 2^ 

X 31 

X 32 

x 33 

X 34 
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The ith row of the matrix h is the row vector £h^ ^ j containing 


all the elements h. . for which i is the first index. 

i, J 


In a like manner the 


jth column of the matrix h is the column vector^ ^containing allthe ele- 


ments h. . for which j is the second index. 

i,J 


The vector indices i = (a,b) and j = (c,d) are equal when their cor- 
responding elements are equal. The addition of indices is also element by element. 


equality of indices: i = j <- 

addition of indices: k = i+j <— 


-> a = c and b = d 
-> k = (a+c, b+d) 


The sum, f, of 2 matrices h and k in this notation is the matrix of 
elements equal to the sum of elements having identical indices. 

f = h + k 


where 




[h. . + k. .] 

i,J 1,0 


The sum of the two pictures x and y is the sample -by-sample sum of 
the pictures . 


Product of Linear Transformations 

The sequential processing of pictures by the relations 
y = xh and then by z = yg 

leads to the product of the matrices h and g since substitution of y into the 
equation for z produces the relation 
z = xhg = xc 


where c = hg. 
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It; terms of elements, the matrix equations are 



_ 

E x n .h. 1 



- "I 

z i,k 


^ y l,j g j,kJ 

substitution 


fz . 


E E Xt ,h. 

l,k 


. . 1,1 1 

J 1 


= 

E . c. 

i 1>1 


where 


[ Ci -k] 


E h 

■j 


i, j S j j 


The product of the matrices is thus a sum over the column index j of the first 
matrix h and. the. row index j of the second matrix g. 

A square matrix h is one which has an equal number of row and column 
indices. The diagonal of a square matrix is the set of elements for which the 
row and column indices are identical. The identity matrix I is then defined to 
he the one which is 1 on the diagonal and 0 elsewhere. 

f 1 !,.:]' [ s i,i] where s 
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Correlation Functions 

The cross -correlation function (D "between the pictures x and z 

xz 

is defined as the matrix product. 

t 

O = x z 
xz 

where the "bar denotes the ensemble average of each element in the matrix product 
x t z. Due to the transposition, there is only one index, 1 = (l,l), over which 
the summation is performed. 


I • 1 


“ ■ — 



! 05 

j xz;;: i> j 


X. Z- 

1, j 


X. z , 

i 4 


It should "be noted that the two pictures x and z can have a different number of 
samples and are not restricted to being the same geometrical size. The cross - 
correlation function of a picture x with itself is known as the auto-correlation 


' 1 

x.x. | 

1 o ! 

L J 

In this 2 -dimensional natation, the correlation functions have many of 
the same properties as the classical one -dimensional correlation functions. The 
transposition of the cross-correlation function between x and z produces the 
cross -correlation between z and x. 


function. 


““ — 


‘ 

<p . . 

XX 1,0 

*— 

X. ,x n . 

i^l 1> <5 


<P. 


xz 


4- 


r t 1 

r — 

1 

t 



) 

X z 

— 


X. ,Z.. . j 

= I Z . , X, . 

1 1,1 Ijj 

u. — 


= Z X = <P 


zx 


The auto-correlation is also its own transpose. 


C p = <D 

^XX XX 
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Statistical independence between the two pictures produces the re- 
sults that the cross -correlation function is the matrix product of the mean 
values of the pictures 


A t -t- 

<p = X Z = X z 

xz 


x and z statistically independent 


where the mean value of the picture x is the row vector made up of the mean 
value of each sample of the picture. In the particular case where the pic- 
tures are statistically independent and one of them has a zero mean, the cross 
correlation function is zero. 


The auto -correlation of the sum of two pictures x and y is the sum 
of the individual auto-correlations and the two cross -correlations. 

(D = (x+y)^ (x+y) = (D+(p+<p+<P 
^ v J J tcx yy xy yx 

Where the pictures are statistically independent and one has a zero mean, the 
auto-correlation of the sum is then just the sum. of the individual auto- 
correlation functions . 

(0 = <p + <p 

xx yy 

The correlation function of a picture made up of independent ele- 
ments all having the same standard deviation a and mean m is an example which 
has many uses. 

2 C 2 

x. x. = cr 8. . + m 
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The correlation function is then 

2 2 

0 = a I + m U 

xx 


where I is the identity matrix and U is the unit matrix of all l's. 


The cross -correlation function between two pictures y^ and y^ linearly 
obtained from pictures x^ and x^ is a matrix operation on the cross -correlation 
functions between and x^» 


" X l h ! + 


*1 


= x 2 h 2+ 


*2 


Vy y = 

^ 1^2 £ 


^ X l h l + V ^ X 2 h 2 + ^2^ 


= hi 


0 

* 1*2 2 


+ h 


t - t 
1 \ 


*2 


+ k.^ x„h 


2 2 


+ k. 




In the case where the pictures x^ and x 2 have zero mean the cross -correlation 

function <p has the particularly simple form 
y l y 2 

0 = ^ „ h 0 + k 1 t k 0 

y^ - 1 2 ^ 

A further simplification occurs in the computation of the auto -correlation func- 
tion 0 of the picture y obtained from the picture x when the picture x is made 

yy 

2 

up of independent, zero-mean elements. In this case, 0 = a I so that the auto 

XX 

correlation function of the picture y is 

0 = o^h + k^k. 

yy 

An example of the use of linear processing theory is in the construc- 
tion of an artificial picture from a random number table. Artificial pictures 
are often used in experiments where it is necessary to know and control the 
properties of the pictures. 
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If it is desired to construct an artificial picture with a prescribed 
auto-correlation function <p , a one-dimensional technique is available for 

yy 

symmetrically factoring the matrix <p to obtain a linear process h which will 

yy 

transform a statistically independent picture into an artificial picture having 
this correlation function. The two-dimensional implications of this technique 
are not presently known. 

A later section describes experiments in artificially creating corre- 
lated pictures from random number tables and linear processing. 


Positive Definite Condition 

An ensemble of pictures x is defined to be linearly dependent whenever 
there exists a nontrivial linear process h of the picture x about its mean x 
which will produce a one-element picture y = x'h having a zero correlation 
function. 



X 1 = X - X. 

A set of pictures in which the. same two samples are always identical 
is an example of an ensemble of linear dependent pictures. When the ensemble of 
pictures is not linearly dependent, it is defined to be linearly independent . 

In the case of an ensemble of linear independent pictures,; every linear process h 
producing a one-element picture has a positive valued correlation function 
t 

= h (fi h > 0 for every h. 

An auto-correlation matrix thus is defined to be positive definite 
whenever it is obtained from an ensemble of linearity independent pictures. 

Where the matrix is positive definite, the determinant \<p | will be non-zero. 

XX 
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The implication of the picture x being linearly dependent is that at 
least one sample x^ is essentially computable from the rest of the picture. By 
definition, if x is linearly dependent, there exists a non -trivial h such that 

y l - S x H h il 


where 


y 


2 


0 . 


This implies that y is essentially zero. 

. ' ; y; .. «• . 0. 


Since at least one of the elements of h has to be non -zero, 
h 5/ 0. This produces the result that 



I 


i^(lA) 



assume that 


Thus, in the case where the auto-correlation matrix has a zero determinant 
|p | = 0, the auto-correlation matrix is non -positive definite so that the 
samples of the picture are linearly dependent . At least one of the samples 
can be produced from the others. In the case of an auto -correlation function 
having a zero determinant |^J = 0, it is always possible to drop the de- 
pendent samples with a linear process h so as to produce a picture w' = x'h with 
a smaller number of linear independent samples and an auto-correlation function 
V w t w i which is positive definite and having a non-zero determinant. 

I^w'w'l ^ °* 

The new picture w' will contain all of the information contained in the ori- 
ginal picture x since the linearly dependent samples are computable from the 
linearly independent samples. 



SPACE -GENERAL CORPORATION 


SGC 203R-U 
Page Ik 


Stationary Correlation Functions 

The initial definition of the correlation function 



permitted each element in the matrix <p to have a different value. In many 

xz 

practical situations the correlation function will he independent of "both trans- 
lation and rotation, and dependent only on the distance of separation between 
the elements x. and z . . 

i 0 

The correlation function <p is defined to be stationary wheh the ele- 
ments of its matrix are functions only of the difference between the indices. A 
stationary correlation function is thus independent of translation but not of ro- 
tation. Whenever <p ^ = £x^z^.j is stationary, there exists P xz (k) such that 

<*> ( k ) = is independent of the index i. The addition and subtraction of 

XZ 1 1 “rK. 

indices is a vector addition and subtraction. A stationary cross -correlation func 
tion satisfies the relation 


M * [vi+k] 


lz A-k 


■ ] 


The stationary auto-correlation function <p satisfies the relation 

XX 
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Symmetric Correlation Functions 


The four elements x^, x^, x^ 2 of the picture 



normally have a correlation function such that 


X 11 

II 

! H 

x 

X 31 

x 32 

horizontal 

X 11 

II 

1 r — 1 

to, 

X 

X 12 

OJ 

KA 

X 

vertical 

X 11 

x 32 


1 ' — 1 
I K\ 

\x 

diagonal 


The first two relations are stationary relations. In the horizontal relation 
the difference in indices (l,2) - (l,l) on the left and {3,2) - (3>l) on the 
right, are the same. The third relation is not a stationary relation since 
the difference of the indices on the left is (2, l) and on the right (2,-1). 

In order to include the diagonal symmetry of the correlation, a sta- 
tionary correlation function is defined to be symmetric whenever the correlation 
function is a function of the magnitude of the elements of the differences be- 
tween the indices. That is, if the stationary correlation function <p(k) is sym- 
metric then <p(k^) = ( p{^) whenever k^ = (a^, b^) and k 2 = (a^, b 2 ) such that 
| a_i j = |a 2 | and |b^| = |b | . When this property holds for the entire correlation 
function, it is said to be symmetric. 
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Distance 

Further restrictions can be placed upon a stationary correlation func- 
tion by specifying that the correlation is only a function of the distance be- 
tween the elements being correlated. The distance d between the picture elements 
having indices i = (i^i^) and j = is defined to be 

<1 = ‘ jp ) 2 + ( i 2 ‘ ^ 2 ) 2 

Two correlation elements are defined to be equal whenever they are the correlation 
function of picture elements separated by the same distance. 

A further modification can be made by specifying the distance measured 
in terms of an ellipse. This modification should have some use in the considera- 
tion of scanned television pictures where the scanning process introduces a dis- 
tortion into the correlation function. 


Stationary and Symmetric Matrix Products 

It would seem intuitive that the matrix product of stationary and sym- 
metric matrices would also be stationary and symmetric. It is fortunate that 
this is the case since It makes it possible to construct invariant symmetric 
processing devices Which are a function -of the distance and direction from the 
particular picture element which is to be. constructed. 


Whenever the matrices h and g are stationary, their matrix product 
c = hg is also stationary. That is if 


c -i 

ik 


? h ij S jk j 

J — 


where h is a function of the difference j-i and g is a function of the difference 
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k-j, then c is a function of the difference k-i. 

c = I S h ■ 2 

i+a,k+a |_ i+a,j j,k+a 

In terms of the differences between the indices k-i = r and k-j 


? ^i4a, j +a,k+a 



s the matrix 


product becomes the familiar convolution formulas. 


c(r) 


E h(r-s) g(s) 


>] 


Sh(s') g(r-s') 

A slightly more useful formular in ccmputation is 

c (r) 


h(r+s ) g(-s)j 
|p(-s) g (r+s )j 


Thus to compute the element c(r), the field of elements of g is reflected through 
the origin and correspondingly multiplied and summed with the field of elements h 
starting at the appropriate element h(r ) . For example where the matrices are 



0 

h 0,i 

0 


"0 

S 0,l 

®1,1 

h = 

h -i,0 

h o,o 

0 

1 — l 

si 

g = • 

0 

S 0,0 

S l,0 


0 

h 0.-l 

0 


0 

w 

0 

0 


the field of g is reflected to give 



“0 

0 

0 ” 

g(-s) 

S l,0 

g 0,0 

0 


s l,l 

g 0,l 

0 


A 
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The term c Q Q of the product is obtained by directly multiplying the 
fields and summing. 


c 


00 


0.0 h Q ± 0 0-0 

h -i,o g i,o h o,o g 0,0 h l, 0 ° 


0 g 


1,1 


h, 


0,-1 s 0,l 


0-0 


h 


0,0 S 0,0 + h -l,0 


"1,0 + h 0,-l S 0,l 


The calculation of the element c of the product is obtained by moving 
the field g(-s) to the element h . , multiplying the fields, and summing. 









°-° h 0,l 

g l,0 

0 g 0,0 

G 11 

= 

h -i,o° h o,o 

g l,l 

h l,0 g 0,l 


i 

°-° V 

H 

o 

0.0 



- 




= 

h 0,0 g l,l + h 0,l 

S l,0 

+ h i,o g 0,l 


Further calculations would produce the matrix 



C 0,2 

%2 


-1,1 

c 0,l 


C 2,l 

- 1,0 

c o,o 

c i,o 

C 2,0 


1— 1 

1 

O 

O 

c i,-i 



where only the non-zero terms have been indicated. 

These are convenient formulas for both hand computation and 
machine computation. 

From the matrix expression it can be shown that the product of two 
matrices h and g which are their own transpose (h^ = h and g^ = g), that the 
product c = hg is not necessarily also its own transpose. 


t 

c 


= g h = gh ^ hg = 


c 


\ 


V \ 
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However, when the matrices are dependent only on the difference between their 
indices and thus stationary, the product will be its own transpose when the 
factors are their own transposes. A stationary matrix h = h(s) is its own trans- 
pose h when h(s) = h(-s). Thus, from the convolution formulas 


c(-r) 


E h(-r-s) g(s) 



s 


= 

E, h(-r+s ' ) g( 


= 

E, h(r-s') g(s 


= 

[c(r)J 


where the first step is the change of variable s ' = -s and the second step 
is the use of the self transpose property on the factors h and g. 

A similar result is that the product of stationary symmetric matrices 
is also symmetric. A stationary matrix is symmetric when it is invariant 
under reflection of any of the elements of the vector representing the index 
difference. That is, where the indices s = (a,b) and s' = (-a,b), the sta- 
tionary matrix g = j^g(s) is symmetric in the first dimension when J^g(s)j = £g(s')j. 
For example 



3 

5 0,0 

2 


0 

1 

0 


is symmetric in the horizontal dimension and not in the vertical dimension. 

Where both the stationary factors of a product have a particular 
symmetry, the product also has that symmetry. In particular where the indices 
are 

s = (a,b) r = (c,d) 

s' = (-a,b) r' = (-c,d) 

the matrix product is 
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where the first step is a change in the order of summation and the second 
step is the use of the symmetry property of the factors h and g. 

It should he noted that fully symmetric stationary matrices have 
the self transpose property. Thus, in the matrix product, the fields can be 
multiplied directly together without a reflection where both of the factors 
are symmetric in both dimensions. 

Construction of a Correlated Picture 

An interesting example of the product of symmetric stationary 
matrices is the processing of a zero-mean picture x whose elements are un- 
correlated to obtain the picture y = xh. The correlation function cp for 

XX 

the picture x is 

2 T 

tn = CT I 

The correlation function <p of the second picture y is then 

yy 

<p = h^ <p h 

yy xx 

= a 2 h t h 

A specific example of a stationary symmetric matrix h is 


w 


w 

1 

w 


w 


Hie matrix multiplication then produces the correlation function 

2 


co = a 

yy 


W 



w 


rp 

< 

ro 

2w 

2w‘ 

2w 

i+4w 2 

2w 

2w 2 

2w 

2w‘ 

. 

2 

w 



w2 


This correlation is then both stationary and symmetric . The 

normalize correlation is plotted in Figure 1 for several values of w. 

o(i) , 

The maximum occurs for w = 1/2. 
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Figure 1, Theoretical Correlation as a Function of 
Distance for various Weighting Coefficients of 
the Linear Process y = xh 


Experimentally the uncorrelated picture of Figure 2 was processed 
with the stationary processor 



1/2 

1 1/2 
1/2 


corresponding to the weighting coefficient w = l/2. This processing produced 
the correlated picture in Figure 3 - The smoothness of Figure 3 is an indica- 
tion of the correlation in contrast to the sharpness of the independent picture 
of Figure 2. 
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Desired Pictures 

Linear processing of the picture x produces another picture y which 
may have one sample or many. Normally, the processing of the picture x is 
directed at obtaining a third picture z, which is known as the desired pic- 
ture, as in Figure k. 


^ . Processed 
Picture y 


Desired y 

Picture z 

Figure Picture Processing Directed at 
Obtaining the Desired Picture z . 

Exact construction of the desired picture z from the picture x is 
normally prohibited by the inherent random differences existing between x and z. 
The error picture in the processing is defined as the difference between the 
pictures y and z . - 

e = y-z. 

The mean squared error of each sample of the error picture is the 
diagonal of the auto -correlation of the error picture 


<P ee = (y-z) t (y-z) 

Optimum Processing 

One approach to the construction of a picture y as close to the pic- 
ture z is to make the squared error between each sample of the two pictures y 
and z as small as possible. In the case where the process is a linear process, 
the process can be represented by the matrix relation 
y = x'h:+:k 


Picture x 
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where h and k are optimum matrices to he specified and x' is the difference 
between the picture x and the mean picture x, x* = x - x. Variation of the 
mean square error 

p ee = (y- z ) 

produces the. .relation 


30 , 


ee 


2 (.ah* x'*-*- dk*) (x’h + k 


- z) 


2dh (V x ^ 9 x , z Y 


+ 2dk (k - z) 


The two relations 


0 , ,h = (ft , 
x x ^x z 


k = z 


are sufficient to minimize each diagonal element in the error correlation matrix 
*P ee ‘ This implies that in the mean square , each element of the produced picture: y 
is.'.. a,s.: close to each element of the desired picture z as can be possible using 
linear processing. • 


The processing which produces a minimum mean squared error for each 
element is thus 


-1 

= x' <p , , <p , + z 

^x'x ’x'z 


This solution assumes that the auto -correlation 0 , , has an inverse. In the 

^x x 

case where the auto-correlation <p , , has no inverse, the determinant \<p , , | 

XX X X 

is zero indicating that it is not positive definite. In this case, there exists 
a linear process h which eliminates a set of samples linearly dependent upon the 
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rest of the samples . The new set of linearly independent samples w contains 
all of the information of the original set of samples . The optimum process 
is then 


y 


w 


CO 1 1 CD | 

r w V r w' 


+ z 


= x' h [hSp , ,h] ^ h^p , 
y x'x' r x'z 


+ z 


where use has "been made of the relation 
w' = x'h 


Simple Example 

A simple example of linear processing results from the consideration 
of pictures z which have had statistically-independent zero-mean noise added 
to them. The picture to he processed is then 

x = z + n 


and the desired picture is z. The auto-correlation of the picture is the sum 
of the auto -correlation functions of the desired picture z and the noise n 




xx 


= 0 + O 

*zz ’nn 


and the cross -correlation <p between the picture x and the desired picture 

xz 

z is 


<p. 


xz 


t 

X z 


t t 
z z + n z 


<P, 


zz 


The mean value of the picture x is equal to the mean value of the desired 
picture z 

x = z . 
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The optimum process is thus, 


y = x' to , . + o J ' 1 «p . . + z 


z z 


nn 


z z 

-1 . - 


where 


= x' (I + , <p )” + z 

' ’z’z r nn 

z’ = z-z and x' = x-x. 


Fpr small values of noise the process is approximately 


y = x - x ' <p~ , , (p 
J z ’nn 


This is a useful result in picture processing since it essentially represents 
a slight "touching up" of the original picture x by the picture -x' <P nn * 

Under normal operation of the transmission facilities, noise is usually small and 
statistically-independent , zero-mean . 

At the other extreme, where the picture is predominantly lost in the 
noise, the optimum processor is 

y = x' o” 1 (p , , + z 
, *nn v z'z' 


In the case of a one sample picture, the processing should be 


y = x' 


'ii 


p u + °n 


+ z 


where 
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In the case of a two sample picture, the process is the following 
operation in terms of classical matrix operations. 




P 22 + a 22 


L'( p !2 + a i2 ^ 


-( p 12 + ° 12 ) 


p ll + °ll' ) 


P 11 P 12 


P 12 P 22 


+ h ± ,z 2 2 


where 


tip (p-q + a i;L ) (p 22 + °22^ " ^ P 12 + CT 12^ 


and 


<p , , = 

^z 'z 


P 11 P 12 


P 12 P 22 


0 H 0 12 
^nn = o 

12 22 


Stationary Stolutiona 

In the case where the experimental picture x and desired picture z 
are stationary, the processing device using the optimum h and k is also a 
stationary device. Under the stationary condition, h is found from the solu- 
tion of the set of equations 


or 


£ ^x'x' ( s_r ) = <Vz (s) 

r 


S ^x'x'^ h ( T+s '> = <P x . z (s) 
r 
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Since there is only one element in the stationary matrix k, the ma- 
trix k is symmetric. Also, since the matrix h is the product of the matrices 
(p , , -1 and <p , , it will he symmetric when the matrices 0 , , and <p , are 
symmetric. In like manner the matrix h will be a function of distance when 
the matrices 0 , , and 0 , are functions of distance. Further simplification 

results in the equations for h where the matrices h, 0 , , and 0 , are funav' 

XX x z 

tions of only distance. All the h's at a particular distance are equal and 
thus require only one equation. That is the set of equations 

z <P x <x.r( s ~ r ) = <P x 'z^ 

needs to be solved for only distinct distances in the index s. 


In particular, for the auto correlation matrix 

Vi V/s 

■Pi f 0 0! 


‘Px'x’ 


and cro&scorrelation matrix 


^x' z 




0), U> 

1 O 1 


a). 


L 

the equations for the linear processor 


h 
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are wh + W,h., = cd 

^00 ni o 


?>l h o + (<P 0 + + 1 = 2 ) h i = “1 


Solution of these equations produce the coefficients 



Generalization of this solution indicates that all coefficients at 
the same distance are equal. The number of simultaneous equations in the solu- 
tion for the coefficients is equal to the number of distinct distances. This 
is a particularly: practical result in that it indicates that the first step 
in the linear processing the picture x is to first add all the samples at a 
particular distance. This produces a computational reduction of about four 
for rectangular grids and up to twelve for hexagonal grids. The hexagonal 
matrix of Figure 5 has a computation reduction of 7 / 2 . 



FIGURE 5 . A Linear Processor Based on a Hexagonal Matrix 
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The coefficients for the hexagonal linear process are 


h = 


^o + 2( ?1 + ^/3 + <° 2 } ~ ^1 


<^>o + ^1 + 2< P/3 + <P 2 ) " 




,axrp - u> o 

1*0 on 


<Po^o + 2{ P± + 2 ^/3 + ^ ■ ^l" 


Minimum Error 


A minimum squared error is produced for' every sample in the linear 


processing 


y = x' h + z 


where 


k ^x'x' ^°x'z 
The use of any other process 


y = x' (h + hg) + (z + kg) 
produces the auto-correlation of the error picture 


<P ee = e e = [x'(h + hg) + (z + kg)-z] [x'(h + hg) + (z + kg)-z] 


where the diagonal entries are the squared sample errors. After some reduction 
with the expression (n , ,h = <p , , the error is 

XX X Zi 
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ee 


= <P ZZ + (z* + kg 1 ) (z + kg) - (£* + k^) z - z 1 (z + kg) 

+ hg 1 ) <p x , x , (h + hg) - (h 1 + hg 1 ) p x , 2 - <p zx ,(h + hg) 

+ (h 1 + hg 1 ) x ' 1 ;(z + kg) + (z 1 + kg 1 ) x' (h + hg) 

= ^z’z' “ ht ^X'x' h + h 5\'x’ h 6 + k 8 \ 


Due to the positive definite condition; the diagonal entries of the term 

hgk<p x , x ,hg and kg^kg are always positive except where hg = 0 and kg = 0. 
Thus, for a linear independent picture x', the mean error is uniquely ob- 
tained by the optimum process 

y = x'h + z 


where 

h = to 1 . o , and x' = x-x. 

^x'x' y x z 

The minimum mean error is the diagonal of the expression 


^Vee = ^z'z' 

. t 

- h '°x 

, ,k 

'x' 

= <p 

z ' z ' 

^zx ' 

^x'x 1 ^x'z 
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In the case of independent additive noise, the minimum error is 


min <p 


ee ^°z ' z 1 <P z 1 z ' ^z ' z 1 + ^nn ^ ^°z ' z ' 


^z ' z ' (<P z ' z ' + ^nn ^ ^nn ' 


-1 


= <p nn (<P z'z’ + ^nn } " <°z’z’ 


= <0 h 
^nn 

The minimum error in the one sample picture of a previous section is 

g 2 _ P 11 C 11 
min P 1;L + P 11 


In the two sample pictures, the minimum error is 

2 

n In n — < 

2 


and 


2 2 

P ll^ g ll g 22" q i2^ + °11 ^ P 11 P 22 ~ P 12 ^ 


'1 min 


2 

a 

'2 min 


(Ppp + g ll^ ^ P 22 + 0 2 ” ^ p i? + 

■ fl io) + °22 ^ P 


P 22^ g ll g 22 “12 


2 *2 

11 P 22 " P 12' 


( p ll + °n ) (P 22 + g 22 ) " (p 12 + g 12 r 


In the case of stationary correlation functions in rotation, the minimum error 
for the two coefficient rectangular grid is 


min <0 = to , , - 03 h - tauh 

^ee *z 1 z ’ 00 11 

and for the two coefficient hexagonal grid, 

min - 03 h - 603 h 

^ee *z ' z ' 00 11 
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Trivial Processing 

The uniqueness of the solution to the optimum processing can be used 
to identify the type of pictures which require no linear processing other than 
possibly a change in scaling. The linear process which is restricted to a change 
in scaling the original picture is 


y = x'd + z 


where d is a diagonal matrix of constants which amplifies each sample independ- 
ently. In the case where the matrix d represents the optimum linear process. 


d = ^x’x’ ^x’z* 

Since d is the unique solution, the only pictures for which a change in scaling 
is the optimum process are where there exists a diagonal matrix d such that 


^®x ' x ' d ^x ' z 

This means that a necessary and sufficient condition is that each column of 
<P X , X , is some multiple of the corresponding column of <P x i z - 

[ith column of <p , ,] d. = [ith column of w , ] every i 
~X XI X z 

The mean error in this case is 
e min = <°z ’ z ’ " d ^’z 


As an example, the processing of a picture z which has added to it correlated 
noise having a correlation function ^ equal to a multiple k of the correla- 
tion function <p , , of the picture z would be processed by a simple change in 
z z 

scaling. In this case, 


aU 


^X’X' 


*^z ' z ' 


+ (D = (l+k) (p , , 
*rm ' ^z'z' 
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and 

*^x ' z ^z ' z ^z ' z ' 
so that the optimum process is 

y = X 'TI*J + 5 

the minimum error is then 

~2 k 

e min 1+k ^z ' z ' 

Orthogonal Pre-processing 

Another approach to linear processing results from the construction 
of an orthogonal maxtrix Q column-wise composed of the normalized eigenvectors • 
of the correlation matrix <p , In particular, 

^x’x’ Q = QX 

where X is the diagonal matrix of eigenvalues X^ corresponding to the eigenvectors 
q . . The normal orthogonal property of the eigenvectors indicates that 

• -T - ~1 T .#?»: V ■ 

Q Q = I or Q = QT 

-Lj near pre-processing of the picture x with the matrix Q produces the 
picture x = xQ which has the auto-correlation function 

X = X q\ - ^ « “ qT Vxx 9 = A 

q q 

and a crosscorrelation function 

t t t _t 
<0 = x z=Qxz=Q (D 

z q *xz 

q 

\ 
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Subsequent optimum linear processing of the new picture x requires that 

q 

, -1 

y = x' X co + z 
J q rx z 

q 


the mean error is then 


2 »-l 

e . = <n , , - (Q A. CO 

min ^z'z' ^zx ^x z . 

q q 

The use of orthogonal processing appears to be of use in pre- 
processing pictures to obtain the picture x^, which is composed of samples 
having zero crosscorrelation and auto-correlation Presently, the experi- 

mental use of this type of pre-processing is limited by the necessity to ob- 
tain the matrix of eigenvectors Q. 

A Stationary Example 

A. particularly useful experimental example of stationary linear pro- 
cessing is the recovery of a correlated picture from a composite picture of 
the correlated picture added to an uncorrelated picture. The composite picture 
could be the result of video transmission of a correlated picture being cor- 
rupted by additive noise, ha a frequency bandwidth - corresponding- to the 
sampling rate of the picture. Another example is the recovery of the uncor- 
related picture from the composite picture of the sum of the correlated picture 
and uncorrelated picture. This second situation corresponds to a high-grade 
picture being corrupted by a fluctuating exposure level resulting in a corre- 
lated bias being added to the picture. 
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The correlated picture to be used in this example is the artificially 
constructed one of a previous 
tionary correlation function 


*1 


= CT 


1/4 


example . The 

picture had 

zero' 

1/2 

1/4 

1 

1/2 


1 

2 

1 l/4 


1/2 

1 

1/2 



1/4 


where w = l/2 was used in order to produce the highest value of normalized cor- 
relation at a distance of 1 sample. The uncorrelated picture to be added to 
the correlated picture also had zero-mean and a stationary correlation function 
composed of only one non-zero element. 



0 

p 

0 a 0 
0 


Since these two pictures are independent of each other, the composite picture 
made from their sum has a correlation function equal to the sum of the corre- 
lation functions and tp . 


<A 


x'x' 


= a 


1/4 



1/4 


1/2 

1 

1/2 

1 

2+a 2 

1 

1/2 

1 

1/2 


1/4 



The composite picture was then used to recover either the correlated picture 
or the uncorrelated picture. 
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Since the pictures are zero-mean and independent, the crosscorrelation 
(P x < z between the composite picture x and the correlated picture is the 

correlation function of the correlated picture. A similar relationship holds 
for the independent picture z^- 

^xz^ ^°i ^xz 2 - ^2 

The matrix equations to be solved are then 

«3 x , x ,h - <p t 1=1,2 

Since the correlation functions , , am are stationary and symmetric > the 

XX 1 "" 

solution h to the equations is also stationary and symmetric. 



By the direct application of the method for matrix multiplication the following 
set of equations is obtained : 


(2 + a 2 ) h Q + 4 h ± + 2 K + h g 

h Q + (13/4 + a 2 ) h x + 2: -h + h 2 

1/2 h Q + 2 h x + (5/2 + a 2 ) h + h 2 

1/4 h Q + h x + h^ + (2 + a 2 ) h 2 



2 

a 

0 

0 

0 
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Only one equation is needed for each of the elements of the matrix h. 
The other equations are identical to these and thus produce no new equations. 
Since the solution of this set of equations is the unique solution to the opti- 
mization of the linear processor, elements of h are zero at greater distance 

than the extent of the correlation functions (p , , and <p , of the experimental 

X X X z 

and desired pictures. In this case 


h . = 0 d > 2 

d 

Experimental Results 

Figure 8 is the composite picture representing the addition of the 

correlated picture, Figure 6, and the uncorrelated picture. Figure J. The addi- 

2 2 
tion was performed with a = 1 so that the correlated picture had a <p oQ = 2a , 

the uncorrelated picture had a <p oQ = a 2 , and the composite had a = 3 cr 2 . 

The correlated picture is thus the dominant picture In the composite picture. 


The recovered correlated picture is Figure 9 and the recovered uncor- 
related picture is Figure 10. As to be expected the dominant correlated picture 
is recovered with better relative precision than the uncorrelated picture. It 
is interesting to note the uncorrelated appearance of the recovered uricorrelated 
picture from the highly, correlated composite picture. 

The minimum recovery error for linear processing of additive pictures 
is the diagonal of the correlation function of the error picture. 

^ee = Pl K + P2 )_1 P2 


= p h 


( 2 ) _ ,( 1 ) 

- <P 2 h 


where h^ 2 ^ is the optimum processor used to recover the second picture and h^^ 
is the optimum processor used to recover the first picture. The recovery error 


Figure 9 - Recovered Correlated Picture y^ from the Composite Picture 


Figure 10. Recovered Uncorrelated Picture y from the Composite Picture 
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is the same for either picture recovered. The correlated and uncorrelated 

2 

pictures were both recovered with a theoretical error of .52a . However, the 

2 

desired correlated picture had a variance of 2a while the desired uncorre- 

2 

lated picture had a variance of a , which indicates a better relative recovery 
of the correlated picture. 


Behavior of the Linear Processor 

A general solution of the set of linear equations of the optimum 
processor can be worked out through some fortunate cancellations and factori- 
zations. For the recovery of the correlated picture, the h's are 


.(1) 


1 


a 2 +l 


+ 


» -? A h (D 


a 2 +l 


h. 


( 1 ) 


a 2 (a^ + 9 /4a 2 + ll/8) 


1 a° + 3l/4a 6 + 131/8 a 4 + 205 /l6 a 2 + 95/32 

h (l) = a 2 (a 2 + 3) (a 2 - 3 A) 

^ 4(a 8 + 31/4 a 6 + 131/8 a 4 + 205/l6 a 2 + 95/32) 


( 1 ) 




a 2 (a 2 + l/2) (a 2 - 3/4) 

6 , /q h Z Ls 2" 


4(a + 31/4 a + 131/8 a + 205 /l6 a + 95/32) 


+ l/2 


a +3 


h 


/2 
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The h's for the recovery of the uncorrelated picture are closely related. They 
are 



2 2^/| 

- 3/4 

2 2 

1 + a a + 1 


h 


1 








As to be expected the sum of the two recovered pictures y^ and y^ is equal to 
the sum x of the correlated picture z^ and the uncorrelated picture z^ . 

y l + y 2 = x ’ h ^ + x ' h ^ 

X,< Px’x'"V + 

= x '<P x i x . 1 (<Pp + <P 2 ) 

= x’ 


= Z 1 + Z 2 


where p x , x , = <£-,_+ <p g 
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The behavior of the h's used to recover the correlated picture is 
plotted in Figure 11. Where the additive corruption is negligible, the pro- 
cessor relies principally on the picture element being processed to produce 
the processed element. (h Q w l, h^ h^> & h^ « 0). At higher levels of 
corruption the processor relies more on the elements near the element being 
processed rather than on the element itself. 

The behavior of the h's in the recovery of the uncorrelated picture 
is plotted in Figure 12. It appears that at all levels of corruption, the 
processor makes use of the elements clustered around the element being processed. 

" ' “2 

For recovery of either picture the theoretical recovery error e is 

computed to be 


e = 


a 2 h^^ 

o 


a* , a 2 (a 2 - 3 A) h (l) 

2 2 n l 

a + 1 a + 1 


2 (l) 

Figure 13 is a plot of this error, ah, normalized by the variance of the 

“ ' ■■ o 

recovered correlated picture. If h is optimum, the variance of the recovered 
picture can be calculated by the formula 


diag<p yy = diag Kz " ^ee 

2 2 

= 2a - a h 

o 

2 2 

The variance of the uncorrelated picture, a a , normalized by the variance of 

2 2 

the composite picture (2 + a ) a is also plotted in Figure 13 as an indication 
of the amount of error present in the composite picture x from which the pro- 
cessed picture y is obtained. 
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FIGURE 1J. Normalized Recovery Error 
Processing the Crater Archimedes 

. Hirther experimental xork was performed on a sampled picture of the 
crater Archimedes, Figure l4. Figure 15 is a composite picture of the crater 
and additive noise of half the variance of the original picture. Optimum lin- 
ear processing then produced Figure 16, which appears to be a decided improve- 
ment over the corrupted picture of the crater. 


I W 


'A' 
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The picture of the crater had the normalized correlation function 




z 


z 1 




.08 


•32 

.54 

•32 

.54 

1 

•54 

•32 

.54 

• 32 


.08 



which is plotted in Figure 17 as a function of distance. 



Figure 17 . Correlation Function of the Crater Archimedes 
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Addition of an independent noise picture with a variance of half that 
of the original picture produces the following set of equations for the deter- 
mination of the optimum processor. 


1.50 h + 2.16 h, 
o 1 

0.54 h + 2.22 h, 

o 1 

0.32 h + 1.08 h, 

0 1 

0.08 h + 0.54 h, 

o 1 

whieh has the solution , 

h = .47 
o 

h-^ = . 12 

Vs - •<* 

hj = -.05 

The large value of h indicates a rather poor suppression of the error. However, 
the experimental results indicate a satisfactory retention of picture features. 
Perhaps, then, mean-square error is not a good indicator of damage done to a 
picture,. by additive noise. 

Suppression of Errors 

The experiment with artificial pictures produced very little suppres- 
sion of independent additive errors. This is to he expected when it is considered 
that the number of processed samples was quite small and the picture was only 
correlated over a small distance. 


+ 1.28 h + O.32 h = 1.00 
/2 e - 

+ 1.08 h^ + 0.54 h 2 = 0.54 
+ 1.66 + 0.64 hg = 0.32 
+ 0.64 h~ + 1.50 h 2 = 0.08 


An upper bound can be obtained for the amount of suppression of addi- 
tive independent noise by assuming a highly correlated desired picture. The 

2 

correlation function of the desired picture z is assumed to be a times the unit 
matrix of all I's. The correlation function of the independent noise is assumed 
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to be a a times the identity matrix of a 1 on the diagonal. All the equations 
for the coefficients of the processing matrix h are identical so that all the 
coefficients are equal. The equation for the single coefficient is then 

(n + a 2 ) h Q = 1 or h = — - — ^ 

n + a 


where there is a total of n coefficients in the matrix h. The variance of the 
error in the optimum recovery is then 


which should be normalized by the variance of the recovered picture 


Thus 




yy 



<°zz " 


ee 


2 

a 


n 


2 

a 

. 2 
n + a 

2 


n + a 


2 


a 


2 

CT 


2 

a 

n 


normalized 


2 2 

In contrast, the variance of the error a a 

2 2 

by the variance (l+a ) a of the picture is 


in the original picture 




nn 


<P. 


xx 


l + a 
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These two relations are plotted in Figure 18 for a processor operating on five 
samples. A suppression of the independent noise by a factor of 5 in variance 
is obtained for small relative levels of noise. As the noise becomes larger 
the linear processing breaks down and produces no improvement in the relative 
quality of the picture. 



Figure l 8 . Theoretical Limits to the Suppression of Additive 
Noise by a Linear Processor Operating on 5 Picture Samples 


The upper bound on the suppression of independent noise indicates 
that the suppression in magnitude of error is inversely proportional to the 
square root of the number of samples. Thus, in order to obtain a suppression 
in magnitude by a factor of 10 the processor requires at least 100 samples to 
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reconstruct one element in the picture. The actual number of samples is pre- 
sently unknown, as the bound given above is computed under the assumption of 
a perfectly correlated picture. Whether or not it is possible to theoretically 
compute the error for a more realistic type of correlation is unknown. The 
computation will require a trick in the inversion of the correlation matrix. 

Minimization with a Constraint 

In the optimum processing of pictures it is often desirable to scale 
the resultant pictures to a particular contrast by having the standard deviation 
of the re sul tant picture equal to the standard deviation of the desired picture. 
Minimization of the mean square error between the resultant picture and the 'de- 
sired picture subject to the constraint that the resultant picture have the same 
standard deviation as the desired picture results in a linear scaling of the 
picture produced under an unconstrained optimization. 

The standard deviation is given by the diagonal terms of the correla- 
tion function of the zero-mean variables. 


diag. 9 , , = diag. <P Z , Z , = constant 


The minimization is then of the diagonal elements of the matrix 


U = (D + \ (0 . . 
r ee v 1 


rr 


where X is an undetermined constant and <0 is the error correlation matrix 

^ee 


P ee = " Z ) (y _ z )' 
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The assumption that the output y is obtained by a linear process on the input 
x requires that 


y = x'h + k y' = x'h 

where x' = x-x and y' = y-y. 

Variation of the matrix U produces the result 

du = s[(y - z^Cy - z) + X y’V' _ 

= 2 j^Bh^ x ,;fe + Bk^)(x'h + k-z) + X . ah* 3 x'^ x'h j 

= 2 »*[*,.*. h(l + X) -(P x , z _ 

+ 2 5k^ j^k — z j 

A sufficient condition for zero variation of the diagonal elements of the matrix 
U is that 

«V>x' h(l+X) 

and k = z 


This is the same solution that was obtained for the unconstrained optimization 
except for a scaling of the output pictures y, producing the desired standard 
deviation. 
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Pre-Processing with Prediction 

One method of transmitting a picture with a sequence of independent 
samples, is to sequentially scan the picture, predict the next element to be 
scanned, measure the prediction error, and transmit the error from which the 
picture can be reconstructed with a similar prediction. 

Theoretically this method will precisely produce a replica of the 
original picture. In those cases where it is possible to do high quality pre- 
diction, considerable reduction can be obtained in the total amount of data 
which must be transmitted . 

Quantization of the transmitted sequence of errors introduces a quan- 
tization drift into the reconstruction process. This drift prevents the system 
from running for any reasonable length of time without having to be reset. 

This section describes a method of prediction and reconstruction from 
quantized errors in such a manner as to avoid the quantization drift in the re- 
construction of the data. 


The picture is to be reconstructed as in Figure 19 where the next 
element in the scan is predicted from the previously reconstructed picture and 
modified by adding the error to this prediction. 


Quantized 

Error 



■> 


Addition 


Recovered 

Picture 



Predictor 




<■ 


Figure 19 . The Reconstruction of a Picture from Quantized 

Prediction Error 
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The prediction error is obtained as in Figure 20 by quantizing the 
error between the actual picture and the predicted picture based upon a recon- 
struction from the quantized error. The picture is first recovered in the same 
manner as it is reconstructed in Figure 19- The prediction of the picture is 
then based upon the reconstructed picture rather than the picture which is to 
be transmitted. The error is obtained by subtracting the predicted picture 
from the actual picture. 



Figure 20, Processing a Picture to Obtain a Quantized 

Prediction Error 

The quantized errors are theoretically a sequence of linearly inde- 
pendent samples. A later section presents a FORTRAN subroutine for transform- 
ing a sequence of independent samples into a binary Huffman Code having a 
maximum entropy per digit. This code essentially minimizes the amount of bi- 
nary bits which must be transmitted in transmitting a picture. Further use of 
various block coding techniques as may be found in Reference 22 will produce 
the necessary reliability in the transmission of the binary data to insure 
reliable reconstruction of the picture. 
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Prediction 

The theory of optimum linear processing provides a means of construc- 
ting a predictor. A preceding Section showed that the optimum linear processor 
for processing the picture x to obtain the prediction of the element z was 

y = x' h + k 

where 

X 1 = X - X 

k = z 

and 

^xlx 1 ? 1 . <°x' z 

A preceding example used the stationary correlation matrix 



The cross correlation matrix between the previously scanned portion of the pic- 
ture and the next sample is then 

1/4 

l/2 l l/2 

1/4 1 X 
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where X marks the center of the matrix. The matrix product to , ,h is obtained 

.X X 

from the formula 


<P x , x . h = E< P X > X .(" r) h(r+s) 

r 

= ^ <P x , x >( r ) h (r+s) 

r 

where the predicting processor has the matrix 



where again X marks the center of the matrix. Equating the matrix product 
to the erosscorrelation produces the following set of equations for the deter- 
mination of the coefficients in the matrix of the processor h. 


+ 

1 — 1 
xl 

OJ 

l/2h 2 

+ 



+ 

S ' 

1 

+ 

1 — 1 
OJ 
1 — 1 

2 h 2 

+ 


+ 

\ 

+ h 6 = 

1 

\ + 

h 2 

+ 

2h 3 

+ 

l/bh k + 

l/2h 5 + l/2hg = 

1/2 


h 2 

+ 

l/** 3 

+ 

2 \ 

1/2K 6 = 

1/2 

h l 


+ 

1/211. 


+ 

2 h 

5 

1/k 


h 2 

+ 

l/2hj 

+ 


+ 2 h 6 = 

1/4 
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Solution of this set of equations would then produce the coefficients for the 
predictor . 


An approximate solution can easily be obtained by assuming that 

h, = 0 so that h, = h_ and h._ = h^. The third order set of equations results 
4 1 2 5 6 H 

in the solution _ _ 


h = 


- 1/2 


- 1/2 

-1/6 2/3 0 

2/3 X 


The power in the prediction error is obtained by the formula 


e 2 = diag (<p z , z , - h* <p x , x , h} 

Quantization of Pictures 

Experimentally a picture is sampled on a regular gird. It is custo- 
mary in high quality pictures to use anywhere from 32 to 128 brightness levels 
in the quantization, of the samples. Due to the limitation in experimental equip- 
ment, the pictures presented in this report are quantized in 6 brightness levels 
(black, 20 $, 40 $, 60$, 80 $, and white). The shade of grey g can be computed 
from the formula 

g = - logyg *> 

where b is the brightness level. The grey scale and brightness levels have the 
relation 


g = 0 1 

2 

3 

4 

5 

6 

7 

8 9 10 

11 

12 

13 

l 4 

, , 1 

1 

1 

1 

1 

1 

1 

111 

l 

1 

1 

1 

b = 1 75 

2 

2/2 
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TE 16/2 32 


U 


128 



SPACE-GENERAL CORPORATION 


SGC 203R-4 
Page 57 


A geometrical pattern is formed by the relative brightness of the 
elements of the picture rather than by the precise brightness levels. The tone 
of the picture is the mean value of the brightness levels. The standard devia- 
tion is a measure of the contrast. In numerous geometrical cases, it is desir- 
able to represent a set of pictures with a uniform mean and standard deviation 
of brightness levels. 

The processing of pictures with linear and non-linear processors pro- 
duces a set of numbers which normally extend beyond the 0-1 range of the bright- 
ness levels-. Several methods are available for quantizing a set of numbers so 
that they can be represented in a finite number of brightness levels . 

One method of representing a picture with quantized brightness levels 
is to distribute the quantization levels so that all levels are represented with 
equal frequency of occurrence. 

Another method of representing a picture is to linearly distribute the 
quantization steps over the range of sample values . This is the method normally 
used in quantizing the samples of an experimental picture. However, pictures 
which are the result of numerical processing, tend to have a few excessively 
small and large values which produce an excessive range. Linearly distributing 
the quantizaticsr steps over this excessive range produces a picture with most 
of ' the quaw-ti yatrl-ffw ' near 50$ • These pictures have very little contrast . 

The quantization of the pictures in this report is based on a linear 
system which distributes the four center shades of brightness (20$, 4-0$, 60$, 
80$) between plus and minus one standard deviation in brightness about the mean 
value of the picture . All values below and above one standard deviation are 



SPACE-GENERAL CORPORATION 


SGC 203R- 4 
Page 58 


made respectively black and white. Figure 21 depicts this type of quantization 
as applied to uniform and Gaussian distributions. For these two distributions 
the mean frequency of occurrence of brightness levels is 

Brightness Level Uniform Distribution Gaussian Distribution 


Black 

.212 

.159 

20 < f > 

.144 

.150 

40# 

.144 

.191 

60# 

,144 

.191 

80$ 

.144 

.191 

Wiite 

.212 

•159 


The frequency of occurrence of each of the six levels is very close to 

l/6 = .167. 

The correlated picture of Figure 3 had the distribution of brightness 
levels shown in Figure 22. Linear quantization between the standard deviations 
produces the results of Figure 3 which clearly show the correlation of the 
picture . 

The principal disadvantage of this method of representation is in the 
comparison of different pictures having different standard deviations. Pictures 
are presented as though they had the same standard deviations. The addition of 
the pictures in Figure 6 and 7 to give that of Figure 8 is an example . " Figure 
6 , 7 , and 8 have relative standard deviations of VF 1 and V 3 • Their standard 
deviations are nomnalized so that a relative comparison of brightness levels can- 
not be directly made. The method, however, does provide a simple way of represen 
ting geometrical patterns with a high degree of contrast in a fairly linear 
manner. 

Figures 23 and 24 are examples of this type of quantization in the rep- 
resentation of the crater Eratosthenes. The white area in the center is an il- 
luminated inner wall which casts the dark shadow extending to the edge of the 
Figure. The measured correlation function is Figure 25 . 
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Figure 21. Quantization of a Picture into Six Brightness Levels in Relation 

to The Gaussian and Uniform Distribution 



Distribution of Real Numbers of the Correlated Picture of Figure 3 
Along with the Quantization Used in its Representation 


Figure 22 . 
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Figure 24. The Crater Eratostenes in 6 Brightness Levels Distributed 
About the Mean Plus and Minus One Standard Deviation 


SPACE-GENERAL CORPORATION 


SGC 203R-4 
Page 6 1 



Figure 25, Stationary Correlation Function of the Crater Eratostenes 


Machine Programming 

The experimental work of this report was performed by hand computa- 
tion. The amount of computation for a very simple picture is at the limit of 
practical hand computation.- 

One of the objectives of this theoretical work has been to cast it 
in a notation which can be easily programmed and run on a large automatic com- 
puting machine. Figure 26 represents the flow charts of the experiment in com- 
puting a linear processor for suppressing undesired errors. Figure 27 is the 
flow chart of a proposed experiment dealing with the removal of redundancy from 
a picture prior to transmission. 
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Figure 27 . Flow Chart of an Experiment in Removing the Redundancy 

from a Picture 
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FORTRAN Subroutines 

The following subroutines have been written to facilitate the experi- 
mental work in picture processing. The FORTRAN convention contained in Refer- 
ence 21 for most IBM 709 and 7090 compilers is used. 

Crosscorrelation of Two Pictures 

The crosscorrelation of the picture X and Y is particularly useful 
in the design of optimum processors. In the following program these pictures 
are assumed to be of the same dimension as in Figure 28. 



X(l,j) and Y(l,j) 

I = 1,M J = 1,N 




Figure 28, The Number of Elements in a Picture 


The correlation to be calculated experimentally 

The dimensions should be odd numbers due to the 

iMD + 1 NC + L 


at the origin, (MO, NO) = - 


has the dimensions of Figure 29 
single value of the correlation 
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C(K,L) 

K = 1,MC L = 1,NC 


Figure 29, The Number of Elements in the Correlation Function 


The crosscorrelation is computed by the formula 
M N 

C(K,L) = jjjj Y / (X(l f J)-XAVG) * (Y(l + K - MO,J + L - NO)-YAVG) 

1=1 J=1 


In those cases where the indices of Y are beyond the range of the picture, the 
average value of Y is used resulting in a zero term being added into the 
summation. 
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FORTRAN Program for Calculating the Crosscorrelation of Two Pictures X and Y 

SUBROUTINE CRCOR (C, MC, NC, X, XAVG, Y, YAVG, M, N) 

DIMENSION X(40,40), Y(40,40), C(40,40) 

DO 1 I = 1,M 

DO 1 J = 1,N 

X(l,J) = X(l,J) - XAVG 

Y(l,J) = Y(l,J) - YAVG 

T = M * N 

MO = (MC+l)/2 

NO = (NC+l)/2 : 

DO 4 K = 1,MC 
DO 4 L = 1,NC 
C(K,L) = 0.0 

DO 3 I = 1,M 
DO 3 J = 1,N 
IY = I+K-MO 
JY = J+L-NO 
IF (IY) 3,3,2 
IF (JY) 3,3,2 
IF (IY-M) 2,2,3 
IF (JY-N) 2,2,3 

2 C(K,L) = C(K,L) + X(I,J) * Y(IY,JY) 

3 CONTINUE 

4 C(K,L) = C(K,L)/T 
RETURN 

END 



Figure JO. Flow Diagram for the Crosscorrelation of 2 Pictures 
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Stationary Matrix Multiplication 
The matrix product 


c = hg 


where c(r) = Jh(s) g(r-s) 
s 

is computed on the assumption that the dimension of each matrix is M x N with 
appropriate subscripts. The center is considered to be the integer part of 


(MO, no) 


M+l N+l ' 

2 ’ 2 


The computation of the product is then 


MH NH 


c(i,j) =y y h(k > l) 


K=1 L=1 


G ( I - MCO + MCO - K + MHO, J - NC0 + NGO - L + NHO) 


where the summation is only over the mutual range of H and G. 
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FORTRAN Program for Calculating the Matrix Product 

SUBROUTINE MATRIX (H, MH, NH, G, MG, NG, C, MC, NC) 
DIMENSION H(40,40), G(40,40), C(40,40) 

MHO = (MS+l)/2 
NHO = (NH+l)/2 
MGO = (MJ+l)/2 
NGO = (NG+l)/2 
MCO =.(MC+l)/2 
NCO = (NC+l)/2 
MO = JCG + MHO - MCO 
NO = NGO + NHO - NCO 
DO 2 I = 1,MC 
DO 2 J = 1,NC 
C(I,J) = 0-0 
DO 2 K = 1,MH 
DO 2 L = 1,NH 
KG = I - K + MO 
LG = J - L + NO 
IF (KG) 2,2,1 
IF (LG) 2,2,1 
IF (KG - MG) 1,1,2 
IF (LG - NG) 1,1,2 

1 C(I,J) = C(I,J) + H(K,L) * G(KG,LG) 

2 CONTINUE 
RETURN 
END 




T 
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Figure 31. Flow Diagram of Stationary Matrix Multiplication 
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CONSTRUCTION OF A HUFEMAN CODE 10 

The conversion of a set of words into another set of words 
having maximum entropy per digit is a coding problem for which an 
exact solution is known in terms of an algorithm. 

The frequency of occurrence of the original words are 
initially listed from the smallest to the largest. The following 
program assumes that the ordering process has already been performed. 

The first two words are distinguished by a 0 for the first and a 1 
for the second. They are then combined and treated as a single word. 

The algorithm proceeds by reordering the new set of words 
distinguishing the first two words with a 0-1 and combining. Each 
word is built up as a sequence of 0 's and l f s depending upon the 
algorithm. 

An index matrix provides the correspondence between the 
ordered set of frequencies and the original set of words. The matrix 

1 10 10 0 

0 0 1 0 0 1 

0 0 0 0 1 0 

indicates that the first, second, and fourth words are associated 
with the first frequency. The third and sixth are associated with the 
second frequency. With this type of index matrix it is then possible 
to know which of the original words are to be distinguished by 
0 's and l*s. 

The program stores the code words reversed from right to 
left with a 1 signifying the end of the word. The code words can be 
obtained from the decimal representation by expressing the representation 
in binary, reversing the ordering and dropping the terminal 1 . 
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A simpler program may be possible by using the concept of "low 
order pairing." A set of low order words are paired and moved up into 
the table without any reordering or moving of the rest of the table. 
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FORTRAN PROGRAM FOR CONSTRUCTING A HUFFMAN CODE 

1 READ 47 N 

2 DIMENSION WORD(N) X(N) FREQ(N) INDEX(N,N) 

3 READ 48 FREQ 

4 DO 10 I = 1,N 

5 DO 10 J = 1,N 

6 IF (i-j) = 9,7,9 

7 INDEX (l,j) = 1 

8 GO TO 10 

9 INDEX (l,j) = 0 

10 CONTINUE 

11 DO 12, I = 1,N 

12 WORD (I) jr.iI.O 

13 L = N 

14 DO 44 M = 2,N 

15 DO 20 J = 1,N 

16 IF (INDEX(1,J))18,18,17 

17 WORD (j) = 2.0 * WORD(j) 

18 IF (INDEX(2,J)) 20,20,19 

19 WORD (j) = 2.0 * WORD(j) + 1.0 

20 CONTINUE 

21 IF (L-2) 45,45,22 

22 DO 25 J = 1,N 

23 X(J) = INDEX(l,j) + INDEX(2, j) 

24 F = IREQ (1) + FREQ (2) 

25 DO 27 K = 3,L 

26 IF (F-FREQ(K)) 27,27,28 

27 CONTINUE 

28 K = K-3 

29 DO 33 I = 1,K 

30 IK = I + 2 

31 FREQ (I) = FREQ (IK) 

32 DO 33 J = 1,N 

33 INDEX (l,j) = INDEX (lK,j) 
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FORTRAN PROGRAM FOR CONSTRUCTING A HUFFMAN CODE (CONTINUED) 

34 K = K+l 

35 FREQ (K) = S 

36 DO 37 J = 1 ,N 

37 INDEX (K,J) = X(J) 

38 K = K+l 

39 L = L-l 

40 DO 44 I = K,L 

41 IK = 1+1 

42 FREQ (I) = FREQ (IK) 

43 DO 44 J = 1 ,N 

44 INDEX (K, J) = INDEX (lK,j) 

45 PRINT 49 (WORD (I)/ I = 1,N) 

46 STOP 

47 FORMAT 

48 FORMAT 

49 FORMAT 
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Figure 52, Flow Diagram for the Construction of a Huffman Code 
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Additional FORTRAN Subroutines 

The following subroutines reirain yet to be written: 

1. Quantization for Display 

This subroutine should compute the quantization constants and 
then quantize the picture, perhaps as it is being read out of the machine. 
Provision should be made for entering the subroutine with the constants already 
given as may be the case in fixed scaling. 

2. Random Picture . Generation 

Several random number routines arc needed for producing artifi- 
cial pictures with known properties. Uniform, and Gaussian distributions with 
zero mean should be the most useful. 

3 . Theoretical Correlation Computation 

This is the computation of the correlation matrix based on the 
linear processing of a picture having a specified correlation matrix. 

4 . Optimum Solution 

Based on the input correlation function <p and the cross- 

x ' x ' 

correlation function <p , , the optimum processor should be obtained. This is 

X z 

the solution of a set of simultaneous linear equations. 
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CONVERGENCE OF ITERATIVE SYSTEMS 

One method of designing optimal processing systems is 
through iterative approximations. An approximation to the system 
is used to compute a Letter approximation iteratively until the 
optimal system is obtained. Sakrison( 2 3) has considered an iterative 
design of a specific optimal system. This section considers some 
of the basic elements of his analysis of the convergence of a system 
to an optimal system. 

NORMS 

In order to consider the convergence of an iterative 
system to an optimal system it is necessary to have some measure of 
the error between the approximation X and the optimum system ©• The 
measure is known as the norm E^. where t denotes the sequence of 
approximations. A very useful norm is the average squared error in 
the case of a single variable or the average inner product of the 
error in the case of a multiple variable system. 

E t = Average £ (X-e).(X-e) ] 

Other norms are based on weighting functions of the error such that 

E t > 0 X f e 

e = o x = e . 


DIFFERENCE EQUATIONS 

The behavior of the norm is of' particular importance in 
describing the convergence of the approximation to optimal systems. 
In those cases where the norm can be described by the first order 
difference equation 

A®t ^ a t + b t E t 



E 


t +1 



quite a bit can be said about the convergence based on the properties 

of the coefficients a. and b . 

x x> 
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The existence of the difference equation which can he 
explicitly solved is the most important idea of this approach. 

Where E* is the solution of the difference equation 

^t = a t + \ E t 

and E. is the actual norm which satisfies the difference relation 

AE , < a . + b . E, 
t — t t t 

the norm E^ is hounded hy the solution E^. 

E t < E t 

and all t whenever 0 < 1 + 

If this were not so, then there would exist a T for which 

Vi > Sfi 

and 

\ ^ 

The difference equations produce the relations 
E T+1 ~ - a T + b T ^ 

^T+l " = ®T + b T E T 

Subtraction produces the inequality 

(Vi • *W s (1+ V (e t • 

which is a contradiction under the assumption that 1 + h^ > 0 since 


Vi -*i+i >0 

E* - E*. < 0 


i 
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Thus for all iterations t, the norm E^ is less than or equal to 
the solution E.* of the difference equation. 

E t < E t 

SOLUTION OF THE FIRST ORDER DIFFERENCE EQUATION 
The first order difference equation 

' a t + \ E t 

has the general solution 

t 


E t+1 


Va. n (1 +b) + E' n Cl + b ) 

^ 1 j=i+l J ° j=o j 


i=o 


In the particular case where hj is the constant b and is the constant a 
raised to the ith power, the solution has the closed form 


E£ - (1+b) 


(gE> 

a 

1 - 1+b 


t+1 


+ Ej, (1 + b) 


t +1 


where a^ = a b^ = b 


Elis solution converges to zero approximately as 


E.J. = ... (1 + b) 


where 


_a_ <1 and 0 < (1 + b) < 1 
1+b 


It is easy to see that if t = 0, the general solution gives the 
correct solution 


E* = a + E‘ (1 + b ) 
1 o o v o' 


SPACE-GENERAL CORPORATION 


SGC 203R-1+ 
Page 80 


* 


By induction, 


the solution for t is assumed correct 


E t 


t-1 

t-1 


n i 

L i 

i=o 

j=i+l 


(1 + V 


t-1 

+ E‘ n (1 + bj 
° j=o J 


From the difference equation the solution for t + 1 is found to be 

E t + i - S + E i t 1 + V 

Substitution of E* produces the relation 


E’ 

t+1 


t-1 



t . 

n (i + O + 

j=i+l J 


E* 

o 


t 

n (1 + b ) 

j=o J 



t 

n 

j=i+i 


(i+ v 


+ 


t 

e* n (i + b ) 
°j=o J 


Which is the relation for the general solution of E^. . 

CONVERGENCE TO ZERO 

Another property of the first order difference equation, 
which makes it useful in the analysis of iterative convergence;, is that 
the solution E* converges to zero whenever the coefficients satisfy 

"C 

the constraints 


0 < 1 + b^ < 1 


oo 



a^ > 0 


OO 


I-i 

=r» 


= A < ® 


(,• ► 
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The logarithm of the second term of the general solution produces 
sequence of relations 


lim In E n (l + b ) = lim 

+,-HCO 0 J=° J "t — ICO 


= lim {in E + Y ln(l + b,)f 

J t- 0 jfe, J 


< In E + lim ) h 


0 & J 


This implies that 


lim E II (l + b.) 
t^co ° J=o J 


Considerably more mathematical rigor is needed to show that the first 
term of the general solution also converges to zero. A sufficient 
proof is one which shows that for every e > 0 there exist a T such that 


V a. n (1 + b.) < e 

L 1 J-i+1 J 


for every t greater than T. 


The proof needs several steps. First, since 0 < 1 + < 1, 

t 

J (l + b ) < 1 and since the a. > 0 there exists some I such that 
+i' y — i - 


a i < 2 


i-T+l 


for every t. 


These inequalities imply that 


l + b j’ - I < 5 


i-I+1 


i-I+1 
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The second part of the proof results from knowing that the product 
converges to zero. That is, there exists some T such that for 
every t > T 


n (1 

j-1+1 


V 


Sa 


This Inequality produces the relation that 


I t 

V a. n (1 + b . ) 

< 

V a n (1 + b ) 
L 1 j=i+i J 

L 1 j«i+i J 

— 

i=o 


i=o 


< !i K 

i=o 


Thus the total sum from 0 to I and I + 1 to t is less than e 
whenever t > T. This implies that the second term converges to zero. 
The rate at which it converges to zero depends on the particular 
coefficients of the difference equation. 


PROPERTIES OF THE ITERATION 


An example of the use of the first order difference equation 
is in the calculation of a norm equal to the inner product of the 
error between the system parameters X and the optimum parameters 8 


E 


n 


Average 



0) • (X. 


n 



The iterative procedure for the determination of the process is 
an iterative correction 


X ,, = X + € 
n+1 n n 
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The error is then 


X 


n +1 


- 6 = X 


n 


9 + 


n 


so that the norm is 


E 


'n+1 


E n + Average 




This is then a difference equation of the norm. 

AE = Average -j e . e + 2(X - 0).c r 

n 0 l n n ' n ' nJ 

In those cases where there exists appropriate a^ and h^ such that 

Average -j e . e + 2(X - 6).e T < a +h E 

°Ln n n ' n) — n nn 

The behavior of the norm E can he analyzed quite effectively with 

n 

the solution of the difference equation. 
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